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OS , ABSTRACT 

We study numerical simulations of satellite galaxy disruption in a potential resembling 
that of the Milky Way. Our goal is to assess whether a merger origin for the stellar halo 
would leave observable fossil structure in the phase-space distribution of nearby stars. 
We show how mixing of disrupted satellites can be quantified using a coarse-grained 
Q\ ' entropy. Although after 10 Gyr few obvious asymmetries remain in the distribution of 

particles in configuration space, strong correlations are still present in velocity space. 
We give a simple analytic description of these effects, based on a linearized treatment 
in action-angle variables, which shows how the kinematic and density structure of the 
debris stream changes with time. By applying this description we find that a single 
dwarf elliptical-like satellite of current luminosity 10 s L disrupted 10 Gyr ago from 
an orbit circulating in the inner halo (mean apocentre ~ 12 kpc) would contribute 
. about ~ 30 kinematically cold streams with internal velocity dispersions below 5 

kms -1 to the local stellar halo. If the whole stellar halo were built by such disrupted 
satellites, it should consist locally of 300 — 500 such streams. Clear detection of all 
these structures would require a sample of a few thousand stars with 3-D velocities 
accurate to better than 5 kms -1 . Even with velocity errors several times worse than 
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O ' this, the expected dumpiness should be quite evident. We apply our formalism to a 

+3 . group of stars detected near the North Galactic Pole, and derive an order of magnitude 

£2 ' estimate for the initial properties of the progenitor system. 
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1 INTRODUCTION 

There have been two different traditional views on the formation history of the Milky Way. The first model was introduced 
by Eggen, Lynden-Bell & Sandage (1962) to explain the kinematics of metal poor halo field stars in the solar neighbourhood. 
According to their view the Galaxy formed in a monolithic way, by the free fall collapse of a relatively uniform, star-forming 
cloud. After the system became rotationally supported, further star formation took place in a metal-enriched disk, thereby 
producing a correlation between kinematics and metallicity: the well-known disk-halo transition. In later studies Searle & Zinn 
(1978) noted the lack of an abundance gradient and a substantial spread in ages in the outer halo globular cluster system. 
This led them to propose an alternative picture in which our Galaxy's stellar halo formed in a more chaotic way through 
merging of several protogalactic clouds. (See Freeman 1987 for a complete review). 

This second model resembles more closely the view of the current cosmological theories of structure formation in the 
Universe. These theories postulate that structure grows through the amplification by the gravitational forces of initially 
small density fluctuations (Peebles 1970; White 1976; Peebles 1980, 1993). In all currently popular versions small objects are 
the first to collapse; they then merge forming progressively larger systems giving rise to the complex structure of galaxies 
and galaxy clusters we observe today. This hierarchical scenario is currently the only well-studied model which places galaxy 
formation in its proper cosmological context (see White 1996 for a comprehensive review). Numerical simulations of large-scale 
structure formation show a remarkable similarity to observational surveys (e.g. Jenkins et al. 1997, and references therein; and 
Efstathiou 1996 for a review). For galaxy formation, the combination of numerical and semi-analytic modelling has proved 
to be very powerful, despite the necessarily schematic representation of a number of processes affecting the formation of a 
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galaxy (Katz 1992; Kauffmann, White & Guiderdoni 1993; Cole et al. 1994; Navarro & White 1994; Steinmetz & Muller 1995; 
Kauffmann 1996; Mo, Mao & White 1998; Somerville & Primack 1999; Steinmetz & Navarro 1999). This general framework, 
where structure forms bottom-up, provides the background for our work. 

We are motivated, however, not only by this theoretical modelling, but also by the increasing number of observations 
which suggest substructure in the halo of the Galaxy (Eggen 1962; Rodgers, Harding & Sadler 1981; Rodgers & Paltoglou 
1984; Ratnatunga & Freeman 1985; Sommer-Larsen & Christensen 1987; Doinidis & Beers 1989; Arnold & Gilmore 1992; 
Preston, Beers & Shectman 1994; Majewski, Munn & Hawley 1994; Majewski, Munn & Hawley 1996). Detections of lumpiness 
in the velocity distribution of halo stars are becoming increasingly convincing, and the recent discovery of the Sagittarius 
dwarf satellite galaxy (Ibata, Gilmore & Irwin 1994) is a dramatic confirmation that accretion and merging continue to affect 
the Galaxy. 

There have been a number of recent studies of the accretion and disruption of satellite galaxies (Quinn, Hernquist & 
Fullagar 1993; Oh, Lin & Aarseth 1995; Johnston, Spergel & Hernquist 1995; Velazquez & White 1995, 1999; Sellwood, Nelson 
& Tremaine 1998). Much of this work has been limited to objects which remain mostly in the outer parts of the Galaxy, which 
may be well represented by a spherical potential plus a small perturbation due to the disk (Johnston, Hernquist & Bolte 
1996; Kroupa 1997; Klessen & Kroupa 1998). In this situation simple analytic descriptions of the disruption process, of the 
properties of the debris, etc. are possible (Johnston 1998). However, it is questionable whether such descriptions can be applied 
to most of the regions probed by past or current surveys of the halo, which are quite local: in this case the influence of the 
disk cannot be disregarded or treated as a small perturbation. 

Since formation models for the Galaxy should address the broader cosmological setting, we are naturally led to ask 
what should be the signatures of the different accretion events that our Galaxy may have suffered through its lifetime. 
Should this merging history be observable in star counts, kinematic or abundance surveys of the Galaxy? How prominent 
should such substructures be? How long do they survive, or equivalently, how well-mixed today are the stars which made up 
these progenitors? What can we say about the properties of the accreted satellites from observations of the present stellar 
distribution? Our own Galaxy has a very important role in constraining galaxy formation models, because we have access to 
6-D information which is available for no other system. Observable structure which could strongly constrain the history of 
the formation of galaxies is just at hand. 

This paper will try to answer some of the questions just posed. We focus on the growth of the stellar halo of the Galaxy 
by disruption of satellite galaxies. We have run numerical simulations of this process, and have studied the properties of 
the debris after many orbits, long after the disruption has taken place. We analyse how the debris phase-mixes by following 
the growth of its entropy and the variations of the volume it fills in coordinate space. We also study the evolution of its 
kinematical properties. In order to model the characteristic properties of the disrupted system, such as its size, density and 
velocity dispersion, we develop a simple analytic prescription based on a linearized Lagrangian treatment of its evolution in 
action-angle variables. We apply our results to derive the observable properties of an accreted halo in the solar neighbourhood. 
We also analyse the clump of halo stars detected near the NGP by Majewski et al. (1994), and obtain an order of magnitude 
estimate for the initial properties of the progenitor system. 

Our paper is organized as follows. Section 2 presents our numerical simulations. In Section 3 we analyse the characteristics 
of the debris in these models, and in Section 4 we develop an analytic formalism to understand their properties. We apply 
this formalism to describe the characteristics of an accreted halo in this same section. In Section 5 we compare our modelling 
with the observations of Majewski et al. (1994). We leave for the last section the discussion of the results, their validity, and 
the potential of our approach for understanding the formation of our Galaxy. 



2 THE SIMULATIONS 

To study the disruption of a satellite galaxy of the Milky Way, we carry out N-body simulations in which the Galaxy is 
represented by a fixed, rigid potential and the satellite by a collection of particles. The self-gravity of the satellite is modelled 
by a monopole term as in White (1983) and Zaritsky & White (1988). 

2.1 Model 

The Galactic potential is represented by two components: a disk described by a Miyamoto-Nagai (1975) potential, 

*n~ ; GMdisk (1) 
v/R 2 + (a + Vz 2 + b 2 ) 2 

where Mdisk = 10 11 Mq, a = 6.5 kpc, b — 0.26 kpc, and a dark halo with a logarithmic potential, 

$haio = ^Liohi(r 2 +d 2 ), (2) 
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Table 1. Orbital parameters for the different experiments. 



Experiment pericentre apocentre 2 max period 
(kpc) (kpc) (kpc) (Gyr) 



1 


10.9 


51.5 


25.0 


0.69 


2 


13.5 


93.1 


69.1 


1.23 


3 


5.0 


51.5 


5.1 


0.64 


4 


9.2 


96.5 


12.0 


1.24 


5 


0.5 


45.5 


30.1 


0.56 


6 


6.0 


37.0 


24.8 


0.48 



with d = 12 kpc and Vhaio = 131.5 km s 1 . This choice of the parameters gives a circular velocity at the solar radius of 
210 kms" 1 , and of 200 kms" 1 at ~ 100 kpc. 

We have taken two different initial phase-space density distributions for our satellites: i) two spherically symmetric 
Gaussian distributions in configuration and velocity space of 1 kpc (5 kpc) width and 5 — 25 kms _1 (20 kms -1 ) velocity 
dispersion, corresponding to masses of ~ 5.9 x 10 7 - 1.5 x 10 9 M Q (4.7 x 10 9 M Q ); and ii) a Plummer profile (1911) 

p(r) = (r2 /°g )5 / 2 - ( 3 ) 

with pq = 3M/47rro, M being the initial mass of the satellite and ro its scale length. In this second case, the distribution 
of initial velocities is generated in a self-consistent way with the density profile. For the characteristic parameters we chose 
M — 10 7 — 10 9 M and r = 0.53 - 3.0 kpc, giving a one-dimensional internal velocity dispersion o\n = 2.9 — 11.3 kms 
The force on particle i due to the self-gravity of the satellite is represented by 

F ( x ') = - (r-* + £ 2)3/2 r " ( 4 ) 

where M ln is the mass of the satellite inside r% = |xi — x c |, x c being the position of the expansion centre defined by a test 
particle with the same orbital properties as those of the satellite. The value for the softening e is 0.25 ro. The approximation 
for the self-gravity of the satellite may not be very accurate during the disruption process, where tidal forces are strong and 
elongations in the bound parts of the satellite are expected. However, because we are interested in what happens after many 
perigalactic passages, well after the satellite has been tidally destroyed, our conclusions on the whole process are unaffected 
by details of the disruption process. 

In total we ran sixteen different simulations, six of which we analyse and describe in full detail in Section 3. Some of the 
remaining simulations are used in Section 4 for comparison with the analytic predictions and the rest are briefly mentioned 
in the discussion. The characteristic properties of our six principal simulations are summarized in Table 1. They differ only 
in their orbital parameters and all initially have a Plummer profile and a mass of 10 7 Mq. We have imposed the restriction 
that the orbits pass close to the solar circle in order to be able to compare the results of the experiments with the known 
properties of the local stellar halo. In all cases the satellite was represented by 10 particles of equal mass. 

In Figure 1 we show projections of orbits 1-6 in three orthogonal planes, where XY always coincides with the plane of 
the Galaxy. Notice that the plane of motion of a test particle on these orbits changes orientation substantially showing that 
the non-sphericity induced by the disk significantly affects the motion of the satellite. 

While orbiting the Galaxy, the satellite loses all of its mass. As expected, the most dramatic effects take place during 
pericentric passages. The satellites do not survive very long, being disrupted completely after 3 passages. This means that for 
our experiments, for any relatively low density satellite on an orbit which plunges deeply into the Galaxy with a period of 1 
Gyr or less, the disruption itself occupies only a relatively small part of the available evolution time. 



3 PROPERTIES OF THE DEBRIS: SIMULATIONS 
3.1 Entropy as a measure of the phase-mixing 

The state of a collisionless system is completely specified by its distribution function /(x, v, t). In making actual measurements, 
it is often more useful to work with the coarse-grained distribution function (/), which is the average of / over small cells 
in phase-space. An interesting property of the coarse-grained distribution function is that it can yield information about the 
degree of mixing of the system (Tremaine, Henon & Lynden-Bell 1986; Binney and Tremaine 1987). 
In statistical mechanics the entropy is defined as 

S = - J d 3 x d 3 v /(x,v,t)ln/(x,v,t). (5) 
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Figure 1. Projections of the orbits of the satellite on different orthogonal planes, where XY coincides with the plane of the Galaxy. All 
distances are in kpc. 



Since the coarse-grained distribution function decreases as the system evolves towards a well-mixed state, an entropy calculated 
using {/) will increase, whereas one calculated using / will remain constant, a consequence of the collisionless Boltzmann 
equation: Df/Dt = 0. We therefore quantify the mixing state of the debris by calculating its coarse-grained entropy as a 
function of time. We represent the coarse-grained distribution function by taking a partition in the 6-dimensional phase-space 
and counting how many particles fall in each 6-D box. Naturally the size chosen for the partition and the discreteness of the 
simulations will affect the result. We can quantify the expected discreteness noise in the following way. The uncertainty in 
the entropy can be attributed to fluctuations in the number counts, which we can estimate as Poissonian, oc \fWl in each 
occupied cell. Therefore, the uncertainty in the entropy in each cell is 



© RAS, MNRAS 000, 



Building up the Stellar Halo of the Galaxy 5 




5 10 0.0 0.2 0.4 0.6 0.8 1.0 1.2 

t (Gyr) 10 t/T m;x 

Figure 2. Evolution of the entropy of the system for the different experiments, as a function of time in (a), and scaled with the mixing 
time-scale in (b). The error in the scaled entropy is of the order of 0.06. 



AS^ — (l + ln-)^ln/V 
for N ^> 1. The total uncertainty is thus 

AS^ (6) 



which, for experiments with 10 s particles is 0.04. In order to have a normalized measure of the mixing properties of the debris, 
we also computed the entropy of points equidistant in time along the corresponding orbit. After a very long integration, the 
orbit will fill the available region in phase-space, whose shape and size are determined by its integrals of motion. In this way, 
by comparing the entropy calculated for the debris with the 'entropy of the orbit', we have a measure of how well mixed the 
debris is. We plot this 'normalized' entropy in Figure 2(a) as a function of time. Note that the orbits which have the shortest 
periods show the most advanced state of mixing, but that this is not complete after a Hubble time. 

The degree of mixing basically depends on the range of orbital frequencies in the satellite, essentially as (Ai/) _1 (Merritt 
1999). This means, for example, that a small satellite will disperse much more slowly than a larger one on the same orbit. 
On the other hand a satellite set close to a resonance will mix on a much longer time scale. One can also imagine that if 
there are fewer isolating integrals than degrees of freedom so that chaos might develop, a satellite located initially in a chaotic 
region will have a large spread Av because of the extreme sensitivity to the initial conditions. Therefore the mixing timescale 
(no longer a phase- mixing timescale) will be very short, since the neighbouring orbits diverge exponentially, instead of like 
power-laws. If indeed the mixing rate is set by the spread in the orbital frequencies v of the satellite, by normalising the time 
variable with this timescale we should be able to derive a unique curve for the entropy evolution S = 5 , max /(t/Tmix). 

In what follows we shall assume that the behaviour of the system is regular as seems to be the case for our experiments. 
Let us recall that any regular motion can be expressed as a Fourier series in three basic frequencies (Binney & Spergel 1984, 
Carpintero & Aguilar 1998). The motion is therefore a linear superposition of waves of the basic frequencies with different 
amplitudes. Terms in this expansion which have the largest amplitude will be the dominant terms and may be used to define 
three independent (basic) frequencies. By performing a spectral dynamics analysis as outlined by Carpintero & Aguilar (1998) 
for ten randomly selected particles in our satellites in each experiment, we compute the frequencies associated with the largest 
amplitude terms in the x- (or y, since the problem is axisymmetric) and z-motions, and their dispersion around the mean. 
We then define 

= m l n{a(^ 1) ),a(^ 2) ),a(^ 1) ),a(^ 2) )}, (7) 
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Figure 3. Isodensity surface of 10 8 po after 14 Gyr, seen from the Galactic plane, for the different experiments. 



The curves obtained by scaling time with T m i x are shown in Figure 2(b) and they can be well fitted with the function 
S = 0.78 - 0.69 exp(-27.03 : 1 



•J max IX 

The good fit and small dispersion confirms that mixing is governed primarily by the spread in frequency. 



(8) 



3.2 Configuration space properties 

To analyse the spatial properties of the debris several Gyr after disruption, we have plotted smoothed isodensity surfaces and 
calculated different characteristic densities. In Figures 3 and 4 we show the density surface at approximately 10 -6 times the 
initial density of the satellite. This encompasses most of the satellite's mass. This density surface practically does not change 
over the last 2 Gyr for experiments 3, 5, 6, showing that the system has reached a stage where it fills most of its available 
3-D coordinate space. The shape of this isodensity surface also gives a measure of how advanced the disruption is. The form 
of the accessible 3-D configuration volume is basically a torus, defined by the apocentre, pericentre and the inclination of the 
orbit. In Figures 3 and 4 we clearly see that shape for experiment 6. Experiments 3 and 5 are in an intermediate state and 
still need to fill part of their tori. In the opposite limit, experiment 2 has filled only a small fraction of its available volume. 
All this is consistent with what was found using the entropy in the previous subsection. The characteristic extent of the debris 
is much larger than the initial size of the satellite. Moreover, debris with these properties may well span a very large solid 
angle on the sky, and so be poorly described as a stream in coordinate space. This is the principal difference between our 
own experiments and those in which the Galaxy is represented by a spherical potential. In the latter the plane of motion of 
the satellite has a fixed orientation, and therefore all the particles have to remain fairly close to this plane, naturally giving a 
stream-like configuration. Late accretion events in the outer halo of the Galaxy will plausibly have this characteristic, as shown 
in Johnston et al. (1996) and Johnston (1998). However, similar behaviour should not be expected in the solar neighbourhood, 
or even as far as 10-15 kpc from the galactic centre since at such radii no strong correlations are left in the spatial distribution 
of satellite particles. Any method which attempts to find moving groups purely by counting stars will probably fail in this 
regime. 
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Figure 4. Isodensity surface of 10 e po after 14 Gyr, seen from the Galactic pole, for the different experiments. 

In Table 2, we present a summary of characteristic densities at different times which were calculated by counting particles 
within spheres of 0.5 kpc radii. The maximum density is achieved at the pericentre of the orbit, though most of the mass is 
distributed closer to the apocentre. In all cases the maximum density is between three and four orders of magnitude lower than 
the initial density of the satellite, and the mean density of the debris is between four and five orders of magnitude lower. These 
values give another estimate of the degree of mixing of the debris. Note that, in accordance with the entropy computation, 
experiment 6 has the smallest characteristic densities, meaning that it has reached a rather evolved state, whereas experiment 
2 has high densities in comparison to the rest. The maximum density in all of the experiments is roughly comparable (similar 
or an order of magnitude lower) to the local density of the Milky Way's stellar halo, though the sizes of regions where this 
density is reached get fairly small, a few kpc 3 , as the evolution proceeds. 



3.3 Velocity space properties 

Let us now focus on the characteristics of the debris in velocity space. We divided the 3-D coordinate space into boxes and 
analysed the kinematical properties of the particles inside each box. Figure 5 shows an example. The scatter diagrams indicate 
that there is a strong correlation between the different components of the velocity vector inside any given box. Notice also the 
large velocity range in each component when close to the Galactic centre. This shows that the debris can appear kinematically 
hot. As we shall see this results from a combination of multiple streams within a given box (clearly visible in Figure 5) and 
of strong gradients along each stream. At a given point on a particular stream the dispersions are usually very small. 



4 PROPERTIES OF THE DEBRIS: ANALYTICAL APPROACH 

In this section we will develop an analytic formalism to understand and describe the spatial and kinematical properties of 
the stream. Let us recall that because the disruption of the satellite occurs very early in its history, the stars that were once 
part of it behave as test particles in a rigid potential for most of the evolution. One of the distinguishing properties of this 
ensemble of particles is that it initially had a very high density in phase-space, and by virtue of Liouville's theorem, this is 
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Table 2. Characteristic densities for the different experiments. 



Experiment time Pmoan Pmax 

Gyr fO 2 M kpc" 3 1O 2 M kpc" 3 



1 5.0 67.0 886.2 

10.0 14.6 223.5 

12.5 7.0 152.8 

15.0 6.8 181.4 



2 5.0 84.7 857.5 

10.0 26.5 376.2 

12.5 11.5 202.4 

15.0 9.7 288.4 



3 5.0 41.5 437.4 

10.0 8.9 72.6 

12.5 8.7 181.4 

15.0 6.9 177.6 



4 5.0 40.8 446.9 

10.0 5.9 99.3 

12.5 5.7 171.9 

15.0 5.1 156.6 



5 5.0 36.4 996.9 

10.0 10.9 210.1 

12.5 6.1 183.3 

15.0 5.7 213.9 



6 5.0 13.8 403.0 

10.0 4.3 82.1 

12.5 4.3 95.5 

15.0 3.4 63.0 



true at all times. At late times, however, this is no longer reflected by a strong concentration in configuration space. This 
evolution can be understood in terms of a mapping from the initial configuration to the final configuration, which we will 
describe by using the adiabatic invariants, namely the actions. 



4.1 Action- Angle variables and Liouville's theorem 

Let H = Zf(q, p) be the (time-independent) Hamiltonian of the problem and (q, p) a set of canonical coordinates. We wish 
to transform the initial set (q, p) to one in which the evolution of the system is simpler, for example, where all the momenta 
Pi are constant. To meet this last condition, it is sufficient to require that the new Hamiltonian be independent of the new 
coordinates H — H(P) — E. The equations of motion then become 

Q t = Ui, P t = 0, 

with solutions 

Qi = Q° + Vit, Pi = P°. 

The generating function that produces this transformation is known as Hamilton's Characteristic function W(q, P), and 
satisfies the Hamilton- Jacobi partial differential equation: 

dW 

The solution to this equation involves N constants of integration a» (including E) for a system with 27V degrees of freedom. 
Therefore, the new momenta P may be chosen as functions of these TV constants of integration. A particularly simple situation 
occurs if the potential is separable in the original coordinate set (q, p). The characteristic function may then be expressed as 
W = J^. Wi(qi, ai...ajv), and the Hamilton- Jacobi equation breaks up into a system of N independent equations of the form: 

„ ( dWi 

tii q%, -= — , cti...a.N 
\ dqi 

each of which involves only one coordinate and the partial derivative of Wi with respect to that coordinate. The transformation 
relations between the original and new sets of variables are 
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Figure 5. Scatter plots of the different velocity components for stars in boxes of ~ 3 kpc on a side at different locations for experiment 
6 at 13.5 Gyr. Similar characteristics are observed in all our experiments. 



P< 



dW _ dW 

dqi ' Ql ~ dPi ' 
and each component of the characteristic function is given by 



Wi{qi,ai..a.N) 



2iPi{q'i,a\..a N )- 



0) 



(For more details, e.g. Goldstein 1953). 

The actions and angle variables are a set of coordinates that describe simply the evolution of a system of particles. They 
are particularly useful in problems where the motion is periodic. The actions are functions of the constants a; and are defined 
for a set of coordinates (q, p) as 

(10) 



(11) 



and their conjugate coordinates, the angles, are 

dW 
<Pi = 7TT- 

OJi 

The evolution of the dynamical system thus becomes: 

fa = 4>° + cii(j)t, 

Ji = J? = constant. 



(12) 



4-1.1 The evolution of the distribution function 

Let us assume that the initial distribution function of the ensemble of particles is a multivariate Gaussian in configuration 
and velocity space 
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/(x, v, t ) = fo exp 



E 



(Xi 



2al 



exp 



E 

i=i 



2a 2 



which we can also express using matrices as 



f(x,v,t ) = /oexp 



1 A t A 

~2 ra a ^^-^ 



(13) 



Here t° denotes the initial time. A^, is a 6-dimensional vector, with three spatial and three velocity components, and A^ 
is obtained by transposing A%. Explicitly A^j = x, — x° for i = 1..3 and A^, 4 = Vj — Vj for i — j + 3 = 4. .6 in a Cartesian 
coordinate system. The matrix a%, is diagonal with cr%u = l/cr 2 for i=1..3, and crj^ = 1/a 2 for i=4..6. As we shall see the 
matrix formulation is particularly useful to study the evolution of the distribution of particles of the system. 

At the initial time, we perform a coordinate change from Cartesian to action-angle variables. Since the particles are 
initially strongly clustered in phase-space, a linearized transformation can be used to obtain the distribution function of the 
whole system in the (0, J) variables. We express this coordinate transformation as 



T°A° 



with 



1 ij 



dw-i 



(14) 



where zu — (x, v), w = (cj>, J) and the elements of matrix T° are evaluated at the central point of the system, around which 
the expansion is performed. By substituting this in Eq. (|l^), and by defining aZ — T° 1 a%T° the distribution function in 
action-angle coordinates becomes 



f(4>, «M°) = /oexp 



1 a t a 
--A„ a w A 



(15) 



that is, it is also a multivariate Gaussian, but with dispersions now given by cr° . 

The deviation of any individual orbit from the mean orbit, defined by the centre of mass or the central particle of the 
system, A Wi — Wi — Wi (t) may in turn be expressed in terms of the initial action-angle variables as 



Ji Ji 

and 



Jj 3 * 



- 4>i(t) = 4>° - 4>° + ^ 1 

aJk 



(Jk - J k )t, 



(16) 



(17) 



where we expanded the difference in the frequencies to first order in Jk — Jk- Eqs. ( |l6[ ) and ( p^ ) can also be written as 
A w (t) = &- 1 (t)A° w , (18) 
where 0(f) is the blockmatrix: 



0(t) 



Is 




-n't 

Is 



(19) 



I3 here is the identity matrix in 3-D, and f2' represents a 3 x 3 matrix whose elements are dQi/dJj. The distribution function 
in action-angle space in the neighbourhood of the central particle at any point of its orbit (4>(t), J) is then 



f(<j>, 3,t) = /oexp \-~Al(t)iT w (t)A w (t) 

with A w (t) = (4> — 4>{t),J — J) and 
a w (t) = ®(tf*l®(t), 

or in terms of the original coordinates a w (t) = (T°0(i)) t a^(T°®(t)). 



(20) 



(21) 



Example: 1-D Case. To understand more clearly what the distribution function in Eq. ( |20[ ) tells us with respect to the 
evolution of the system, we consider the 1-D case. The initial distribution function becomes: 

^-^ f (J — J) 2 



f(<f>,J,t ) = /oexp 



2-5 



<t?)(J - J)C 0J 



where C^j denotes the initial correlation^] between <f> and J. After considering the time evolution of the system (as in Eq. (|l7 
we find 



r' ~2 'l 

Cj,j is not the correlation coefficient, usually denoted as p. They are related through p = — — - — •> \> . 

1 - c 4>J rT 4,' T j 
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Figure 6. 1-D graphical interpretation of Liouvillc's theorem and the evolution of the system in phase space. The system is initially 
a Gaussian in action-angle space, with no correlations between tf> and J. As time passes by, the system evolves into an ellipsoidal 
configuration, with principal axes that are no longer aligned with the action or the angle directions. After a some time, the system wraps 
around in the angles, giving rise to phase-mixing: at the same phase we observe more than one stream, each with a small variance in the 
action due to the conservation of the area in phase-space. 



f(<j>, J, t) = f exp 



(J- 



2c? 



4>{t)){j-j) 



W,J + — 2" 



where Q' — dQ,/dJ. This means that the dispersion in the J-direction effectively decreases in time and the covariance between 
<j> and J increases with time. The system becomes an elongated ellipsoid in phase-space as time passes by as a consequence of 
the conservation of the local phase-space density. This evolution is illustrated in Figure 6. 



4-1.2 The distribution function in observable coordinates 



To compute the characteristic scales of a system that evolved from an initial clumpy configuration, such as satellite debris, we 
have to relate the dispersions in action-angle variables to dispersions in a set of observable coordinates. The transformation 
from the action-angle coordinate system to the observable (x, v) has to be performed locally since we generally cannot express 
in a simple way the global relation between the two sets of variables. Because the system has expanded so much along 
some directions in phase-space, the transformation from (</>, J) to (x,v) has to be done point to point along the orbit. This 
transformation is given by the inverse of T at time t: 



(22) 



where the derivatives are now evaluated at the particular point of the orbit around which we wish to describe the system in 
(x,v) coordinates. In particular, if the expansion is performed around ((j>(t),J) then 



A„(t) 



>(*). 



and the distribution function may be expressed in the region around vj = (x, v) as 
/(x,v,i) = /oexp -iA ro (t) t a ro (t)A ro (t) 
with 

A roi (t) = 
and 

a m {t) = (T°0(t)T- 1 ) t ( £(T o 0(i)T- 1 ) 



Xi(t), 



i = 1..3, 
j +3 = 4.6, 



(23) 



(24) 



(25) 



(26) 
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We find once more that, locally, the distribution function is a multivariate Gaussian, where the variances and covariances 
depend on their initial values, on the time evolution of the system and on the position along the orbit where the system centre 
is located at time t. 

If we wish to describe the properties of a group of particles that are located at a different point w than the central particle 
(i.e. the expansion centre does not coincide with the satellite centre at time t) a slightly different approach must be followed. 
The region of interest is then A w (i) = w' — w(t) = (w' — w) + (—w(t) + w) = A^ + D(t). We replace this in Eq. ( pc| ) and 
write 

1 

2 



f(<p,J,t) = /oexp 
or equivalently 
f((f>,J,i) = /o(t)exp 



{A' w - t)(t)y a w [t) (A' w -T)(t)) 



(27) 



■iA" 



a w (t)A' w - ^A'J a w {t)t)(t) 



(28) 



where fb(t) = fo exp [— 1/2 T)(t)^cr w (t)T)(t)]. We may now express A^, = T' 1 A^ 7 , since the transformation is local again. 
The distribution function becomes 
1 



/(x',v',t) = /o(*)exp 
with 

5(t)=T'D(t), 



|(A^ -«(*)) W(t)(A^-«(t)) 



a^(t) = (T'- 1 ) t a u ,(t)T' 



(29) 



(30) 



and fo(t) — fo(t) exp [— 1/2 (T~ 1 S(t)ya w (t)T~ 1 S(t)]. This means that the local distribution function is Gaussian centered 
around x m = x + <5(t), which in general will not be very different from x, with variances given by the elements of cr ro /(t). Thus 
the same type of behaviour as derived for the region around the system centre holds also if far from it. 

The formalism here developed is completely general, but the actions will not always be easy to compute. As we mentioned 
briefly in the beginning of this section, this depends mainly on whether the potential is separable in some set of coordinates. 
We focus on the spherical case and a simple axisymmetric potential in the next section to show how this procedure can be 
used to describe the characteristic scales of the debris. We refer the reader to the Appendix for details of the computation. 



4.2 Spherical Potential 

4-2.1 Analytic predictions 

For a spherical potential 3>(r), the Hamiltonian is separable in spherical coordinates and depends on the actions J v and Je 
only through the combination J v + Je = L. This means that the problem can be reduced to 2-D, and so we may choose a 
system of coordinates which coincides with the plane of motion of the satellite centre. The position of a particle is given by 
its angular (ip) and radial (r) coordinates on that plane. Thus 

L = J.0 = Pi/,, 

J r = - [ dr-J2{E-${r))r 2 - L 2 , (31) 
7r / r 

where L is the total angular momentum of the particle, E its energy and ri and ri are the turning points in the radial 
direction of motion. The frequencies of motion and their derivatives needed to compute the matrix &(t) and to obtain the 
time evolution of the distribution function, can be obtained by differentiating the implicit function g = g(E, L, J r ) = defined 
by Eq. d§). 

Let us assume that the variance matrix ^] in action-angle variables is diagonal at t = 0. This simplifies the algebraic 
computations and, since we are only trying to calculate late-time behaviour, this assumption does not have a major influence 
on our results. As shown in the previous section, the evolution of the system in action-angles is obtained through a w (t) — 
&(t) *a° w &(t). We find the properties of the debris in configuration and velocity space by transforming the action-angle 
coordinates w = (cj>, J) locally to the separable ui = (x, p), and then by transforming from uo = (x, p) to m = (x, v). That is 
o- ro (t) = T' t cr I1 ,(t)T', with the T' = T m ^T p ^ v . 

The diagonalization of the variance matrix a m (t) yields the values of the dispersions along the principal axes and their 
orientation. It can be shown that two of the eigenvalues increase with time, whereas the other two decrease with time. This is 
directly related to what happens in action-angle variables: as we have shown for the 1-D case, the system becomes considerably 
elongated along an axis which, after a very long time, is parallel to the angle direction. For 2-D (3-D), the evolution in action- 
angles can also be divided into two (three) independent motions (whether or not the Hamiltonian is separable) , so that along 



t Strictly speaking a is the inverse of the covariance matrix. However we will loosely refer to a as the variance matrix. 
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each of these directions this same effect can be observed. The directions of expansion and contraction are linear combinations 
of the four axes (e^, e r , , e Vr ) and, generally, none is purely spatial or a pure velocity direction. 

To understand the properties of the debris in observable coordinates, we will examine what happens around a particular 
point in configuration space. This is equivalent to studying the velocity part of the variance matrix: cr ro (w). For example, by 
diagonalising the matrix cr ro («) we obtain the principal axes of the velocity ellipsoid at the point x. Its eigenvalues are the 
roots of det[a-^(v) — AX] = 0. For t ^> t or b 



A1A2 
Ai + A 2 



+ 4 in' O' O' 2l \ 2 „, 2 P r 

t ("ll"33 - "13 ) r 7^20-11033, 



0"11 



fill 



^3 



+ f 



f 2 1 



+ 033 ni 



^33 



1/a 2 ,., the initial variance in the angles. Since 



where the subindices 1 and 3 represent if) and r respectively, and an 
a(vi) = yl/Ai both directions in velocity space have decreasing dispersions on the average. 

So far we did not describe how the debris is spread along the transverse direction to the plane of motion: and e v # . This 
is because we reduced the problem to 2-D in configuration space. However, the problem is not really 2-dimensional since the 
system has a finite width in the direction transverse to the plane of motion. Now that we have understood the dynamics of 
the reduced problem, the generalization to 3-D is straightforward. If the variance matrix initially is diagonal in action-angle 
variables, then the dispersions along <f)$ and J$ do not change because the frequency of motion in the transverse direction 
is zero. Thus the velocity dispersion and width of the stream also remain unchanged in the direction perpendicular to the 
orbital plane. 

By integrating Eq. (E4) with respect to the velocities, we compute the density at the point x 



p(x,t) 

For t > iorb, 



dvg dv v dv r /(x, v, t) . 



p(x,t) 



(27r) 3 / 2 /oo-^ 

1^11^33 ~~ ^13 




r 2 sin #|p r pg| t 2 ' 



(32) 



(33) 



where ax is the initial dispersion in the quantity X. This equation shows that the density at the central point of the system 
decreases, on the average, as 1/t 2 . It tends to be larger near pericentre since it depends on radius as 1/r 2 ; moreover it diverges 
at the turning points of the orbit. Even though the system evolves smoothly in action-angle variables, when this behaviour 
is projected onto observable space, singularities arise associated with the coordinate transformation. In action-angle variables 
the motion is unbounded, whereas in configuration space the particle finds itself at a 'wall' near the turning points. This 
divergence shows up in the elements of the transformation matrix T w _nn (Eq. (A3)), some of which tend to zero, while others 



diverge keeping the matrix non-singular. Because of the secular evolution of the dispersions, the intensity of the spikes will 
decrease with time. They are generally stronger at the pericentre of the orbit than at the apocentre, because of the 1/r 2 
dependence of the density. 

A direct consequence of the secular evolution is that the characteristic sizes of the system, the width and length of 
the stream, will increase linearly with time, reflecting the conservation of the full 6-D phase-space density. At the turning 
points one of these scales becomes extremely small. In Figure 7 we plot the predicted behaviour of the dispersions along the 
principal axes of the velocity ellipsoid as a function of time. We have chosen for the initial conditions a spherically symmetric 
Gaussian in configuration and velocity space. We follow the evolution of the variance matrix and, in particular, of the velocity 
dispersions along the three principal axes at the positions of the central particle. In all panels we can clearly see the periodic 
behaviour associated with the orbital phase of the central particle, superposed on the secular behaviour related to the general 
expansion of the system along the two directions in the orbital plane. The dispersion in the third panel is on average constant: 
it is in the direction perpendicular to the plane of motion. Its periodic behaviour is due to the fact that we did not start 
with a diagonal matrix in action-angles. The initial transformation from (x, v) to ((f), J) produces cross terms between all 
three directions. As the system evolves, and we project again onto configuration space, our 6-D ellipsoid rotates continually, 
producing a contribution in the direction perpendicular to the orbital plane which varies with the frequencies fl r and Qe- By 
fitting a(v) /ao(v) = a/(l + t/to), we find for the velocity dispersion in the first panel a = 1.5 and to = 0.6 Gyr, whereas for 
the dispersion in the second panel a — 2.6 and to —0.1 Gyr. 

In th e last panel we show the behaviour of the product of the three dispersions, which is proportional to the density (see 
Eq. (A1C)). Note that, since two of the velocity dispersions have decreased approximately a factor of ten, the density has done 
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Figure 7. Time evolution of the velocity dispersions along the major axis, computed as outlined in Section 4.2, for the logarithmic 
spherical potential of Eq. (|^). Two of the dispersions decrease with time as 1/(1 +t/to) (dotted curve), whereas the third one is constant 
on the average. The periodic variations are due to the combination of the radial and angular oscillations, as described in the text. The 
last panel shows the product of the three dispersions which is proportional to the density (full curve). The radial oscillation is shown 
(dotted curve) so that the occurrence of density spikes can be compared with the location of the turning points of the orbit. 



so by a factor of hundred. Note also the decrease in the amplitude of the spikes and the good correlation of these with the 
turning points of the orbit. 



4-2.2 Comparison to the simulations 

In order to assess the limitations of our approach, we will compare our predictions with simulations of satellites with and 
without self-gravity. We first consider what happens to a satellite with no self-gravity moving in a spherical logarithmic 
potential. We take two different sets of initial properties for the satellites: lkpc width and o\d = 5kms _1 , corresponding to 
an initial mass of ~ 5.9 x 10 7 Mq; and 5kpc width and am = 20kms~ 1 , corresponding to M ~ 4.7 x 10 9 Mq for the larger 
satellite. Both begin as spherically symmetric Gaussians in coordinate and velocity space. We launch them on the same orbit 
so that we can directly study the effects of the change in size. 

What observers measure are not the velocity dispersions or densities of a stream at a particular point, but mean values 
given by a set of stars in a finite region. We can estimate the effects of this smoothing by comparing our analytic predictions 
with the simulations. In the upper panel of Figure 8 we show the time evolution of the density (normalized to its initial value) 
for the small satellite. The full line represents our prediction and the stars correspond to the simulation. We simply follow 
the central particle of the system as a function of time, and count the number of particles contained in a cube of 1 kpc on a 
side surrounding it. Triangles represent the number density from an 8 times larger volume (2 kpc on a side). The agreement 
between the predictions and the estimated values from the simulations is very good. The representation of a continuous field 
with a finite number of particles introduces some noise which, together with the smoothing, is responsible for the disagreement. 
Note, however, how well the simulated density spikes agree with those predicted at the orbital turning points. The overall 
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Figure 8. Time evolution of the density for a satellite moving in a spherical potential (Eq. (^)), with similar orbital parameters as 
those of Experiment 6 in Table 1. The full line represents our prediction, normalized to the initial density. In the upper panel we plot 
the density behaviour for the ~ 5.9 X 10 7 Mq satellite (see main text), whereas the lower panel corresponds to the ~ 4.7 X 10 9 Mq 
satellite. The stars indicate the number of particles that fall in a volume of 1 kpc on a side around the central particle of the system, and 
the triangles represent the number of particles in a cubic volume of twice the side, both normalized to the initial value. The spike-like 
behaviour occurs at the turning points of the orbit (see main text — Eq. (B3ft). 



agreement is slightly better for the small cube than for the large one. This is due to the smoothing which inflates some of the 
dispersions as a result of velocity gradients along the stream. 

In the lower panel of Figure 8 we show a similar comparison for the large satellite. In general the prediction does very 
well here also. Note for the small boxes and at late times, we only have simulation points at the spikes (i.e. when the density 
is strongly enhanced). This is because the satellite initially has a larger velocity dispersion and therefore spreads out more 
rapidly along its orbit. 

We tested the effect of including self-gravity in the small satellite simulation, and found no significant qualitative or 
quantitative difference in the behaviour. 



4.3 Axisymmetric case 

As an illustrative example of the main characteristics of the axisymmetric problem, let us consider the class of Eddington 
potentials $(r, 0) — $i(r) + r]((3 cos 9)/r 2 (Lynden-Bell 1962; 1994) which are separable in spherical coordinates. The third 
integral for this type of potentials is Iz — ^L 2 + i](f3cos8). The actions are computed from: 

J v = L„ (34) 



J e = n= HV 2 ( f >-lW)-773. ( 35 ) 
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Figure 9. Time evolution of the velocity dispersions along the principal axes, computed as outlined in Section 4.2 and 4.3, for the simple 
axisymmetric potential of Eq. (^). Now all the dispersions decrease with time as 1/t (dotted curve). The periodic time behaviour is 
due to the combination of the radial and angular oscillations, as described in the text. The last panel shows the product of the three 
dispersions which is proportional to the density. The radial and S-oscillations are also plotted to indicate the position of the turning 
points. 

Since the frequencies of motion are all different and non-zero, the system has the freedom to spread along three directions 
in phase-space. The conservation of the local phase-space density will force the dispersions along the remaining three directions 
to decrease in time. 

Following a similar analysis as for the spherical case we derive for the density at the central point x(t) of the system at 
time t 

(2^) 3/2 /o 



p(x,t) = 



dh 



ydet^o |det«'| dJ e r 2 sin 6\p rPg \ t 3 : 



(37) 



where is the angle submatrix of the initial variance matrix in action-angle variables. Therefore the density at the central 
point of the system decreases as t~ 3 , because of the extra degree of freedom that the rupture of the spherical symmetry 
introduces (see Appendix B), and so after a Hubble time, the density decreases by approximately a factor of a thousand. 

In Figure 9 we plot the time evolution of the components of the velocity ellipsoid for a system on an orbit with the same 
initial conditions as for the spherical case, in the potential 



$(r,0) = ^log (r 2 



/3 2 



(38) 



where «h = 123 kms - , d = 12 kpc and /3 = 950 kpc kms . This choice of parameters produces a reasonably fiat potential 
which is physical (giving a positive density field) outside 7 kpc. All velocity dispersions now decrease as 1/t. 

The analytic formalism d evel oped here can be applied to any separable potential in a straightforward manner, using the 
definitions and results of Sec. 4.1. This includes, of course, the set of Stackel potentials which may be useful in representing 
the Milky Way (Batsleer & Dejonghe 1994), or any axisymmetric elliptical galaxy (de Zeeuw 1985, Dejonghe & de Zeeuw 
1987). The only difference is that the matrix T of the transformation from the usual coordinates (x,v) to the action-angle 
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variables should be first multiplied by the matrix of the mapping from (x,v) to the ellipsoidal coordinates (A, jj,, ip,p\,p^,p v ), 
since this is the system in which the problem is separable. We discuss some of the properties Stackel potentials and derive, 
for a particular model for our Galaxy, the explicit form for the density in Appendix C. Even if the potential is not separable 
our general results on the evolution of the system remain valid provided most orbits remain regular. In the general case the 
frequencies and their derivatives with respect to the actions will have to be computed through a spectral dynamics analysis 
similar to that used in Section 3.1 (Carpintero and Aguilar 1998). 



4.4 What happens if there is phase-mixing 

The procedure outlined above assumes that only one stream of debris from the satellite is present in any volume which is 
analysed. When phase-mixing becomes important we may find more than one kinematically cold stream near a given point. 
The velocity dispersions of the debris in such a region would then appear much larger than predicted naively using our 
formalism. We can make a rough estimate for the velocity dispersions also in this case by using the following simple argument. 

If the system is (close to) completely phase-mixed, then the coarse-grained distribution function that describes it will be 
uniform in the angles and therefore will only depend on the adiabatic invariants, i.e. /(x, v) = /(J(x, v)). Since these are 
conserved the moments of the coarse-grained distribution function will be given by the moments of the initial distribution 
function. Therefore /(J) is completely determined by the initial properties of the system in the adiabatic invariants space. If 
the initial distribution function is Gaussian in action-angles then /(J) will be Gaussian with mean and dispersion given by 
their values at t — t°. 

As an example, let us analyse the velocity dispersion in the (^-direction in a particular region in which there is a multistream 
structure: 

r , , , ,9 , , u f d 3 x d 3 J (2*. - £t) /(J) 
2( v = Jd 3 xd 3 v (vy-vyf /(J(x,v)) = I \R R) M 

° [Vv) fd 3 x d 3 v /(J(x,v)) Jd 3 xd 3 Jf(J) 

where we used that v v = J v /R. By expanding to first order we find 

a 2 {v^) = o 2 {J v )/R 2 + A 2 X JI/R\ (39) 

Here we replaced <j(R) by (the size of the region in question) which is justified by our previous result that the spatial 
dimensions of streams grow with time; and neglected the correlation between J v and R. The first term in Eq. (^) estimates 
the dispersion between streams, while the second estimates the contribution from the velocity gradient along an individual 
stream. For the experiments of Table 1 the values of the dispersions range from 50 to 150 kms -1 . These dispersions increase 
in proportion to those of the initial satellite. 



4-4-1 The filling factor 



We can use the results of our previous section to quantify the probability of finding more than one stream at a given position 
in space. This probability is measured by the filling factor. We define this by comparing the mass- weighted spatial density of 
individual streams with a mean density estimated by dividing the mass of the satellite by the total volume occupied by its 
orbit. The first density can be calculated formally through an integral over the initial satellite: 



{p{t)) = ~M dm ( x ' v ) P( x » v )(*) = Jf 



d 3 x d 3 vf{x,v,t°) p(x,v)(t), 



where p(x,v)(t) is the density at time t of the individual stream in the neighbourhood of the particle which was initially at 
(x, v) . The filling factor is then 

M 1 



F(t) 



Vo{ P {t)Y 



where V is the volume filled by the satellite's orbit. An estimate of the filling factor can be obtained by approximating {pit) ) 
by / 9(x,t)/(2 v / 2) taken from Eqs. @, (0) or ( |ci^ ) for spherical, axisymmetric Eddington or Stackel potentials respectively. 
The factor \j1\pl is the ratio of the central to mass- weighted mean density for a Gaussian satellite. We approximate V = 
4?r r 3 po cos 8{ /3, where r apo and 8{ correspond to the orbit of the satellite centre. Since we are interested in deriving an estimate 
for the filling factor for the solar neighbourhood, we focus on the Stackel potential described in Appendix C, which produces 
a flat rotation curve resembling that of the Milky Way. Thus 



F(t) 



6V2M v /det^> <#)< 



Idet n'l 



2(2tt)5/2 f 



o cos 9i 



dh 



dh 



(40) 
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where A, v are spheroidal coordinates (for which the potential is separable), J\ and J„ are the corresponding actions, and 



Q,\ and f2„ the frequencies; and 73 is the third integral of motion. If we approximate (v\V v ) 
M/(2na(x)a(v)) 3 then 



„/4 and replace /o 



where 



Corbit — 



3y/n (\u-X\) (R) v* lTC |detn'| 



2 cos #f 



TTT 



dJx dJ v 

depends on the orbital parameters of the satellite, and 



Cic 



with 



A 



R 



ft T = 2p T 



dh 
dJ x 

dp T 
Or ' 



dh 
8J V 



P S Q 3 



(41) 



(42) 



(43) 



A, v, 



is a function of its initial position on the orbit. (See Appendix C for further details and definitions). This last expression holds 
if the satellite is initially close to a turning point of its orbit. 

For example, a satellite of 10 km s -1 velocity dispersion and 0.4 kpc size on an orbit with an apocentric distance of 13 
kpc, a maximum height above the plane of 5 kpc and an orbital period of ~ 0.2 Gyr, gives an average of 0.4 streams of stars 
at each point in the inner halo after 10 Gyr. A satellite of 25 km s _1 dispersion and 1 kpc size on the same orbit would produce 
5.9 streams on the average after the same time. Let us compare this last prediction with a simulation for the same satellite 
and the same initial conditions in the Galactic potential described in Section 2. In Figure 10 we plot the behaviour of the 
filling factor from the simulation, computed as 

where N is the total number of particles, nit) = N^ 1 £\ pi with pi the density of the stream where particle i is, which we 
calculate by dividing space up into 2 kpc boxes and counting the number of particles of each stream in each box. Note that 
the filling factor increases as t 3 at late times as we expect for any axisymmetric potential. Our prediction is in good agreement 
with the simulations, showing also that it is robust against small changes in the form of the Galactic potential. 



4.-4-2 Properties of an accreted halo in the solar neighbourhood 

To compare with the stellar halo it is more useful to derive the dependence of the filling factor on the initial luminosity 
of a satellite. We shall assume that the progenitor satellites are similar to present-day dwarf ellipticals, and satisfy both a 
Faber- Jackson relation: 

log A _ 3.53 log -p^- ~ 2.35, (44) 
Lq kms" 1 

for Hq ~ 50 km s _1 Mpc _1 , and a scaling relation between the effective radius (R e ~ &(x)) and the velocity dispersion o~(v): 

log j^pr- 1-15 log jA- 1-64, (45) 

both as given by Guzman, Lucey & Bower (1993) for the Coma cluster. Expressed in terms of the luminosity of the progenitor, 
the filling factor then becomes 

f T \ °' 776 

F(t) ~ C orbit Cic \j-J (fl A t) 3 , (46) 
where L n is a normalization constant that depends on the orbit and on the properties of the parent galaxy as: 

rapo ^ ( ^ V' 29 . (47) 

10 kpc/ V 200 kms- 1 / v ; 



L n = 3.75 x 10 n L Q 



If the whole stellar halo had been built from disrupted satellites, we can derive the number of streams expected in 
the solar neighbourhood by adding their filling factors using the appropriate orbital parameters in Eq. @ or Eq. @: 
i*o(i) = N sa .tF(t). For a sample of giant stars located within 1 kpc from the Sun with photometric distances and radial 
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Figure 10. Time evolution of the filling factor for a satellite with an initial velocity dispersion of 25 kms - 1 and size of 1 kpc, moving 
in the Galactic potential described in Section 2. Its orbital parameters resemble those of halo stars in the solar neighbourhood. The 
dashed-curve indicates a 70 + 71 f 3 fit for late times. 



velocities measured from the ground (Carney & Latham 1986; Beers & Sommer-Larsen 1995; Chiba & Yoshii 1998), and 
proper motions measured by HIPPARCOS, we estimate C or bit X Cic ~ 1-29 x 10 -3 . The median pericentric (apocentric) 
distance is 3.7 (11.6) kpc, and the median Q,\ is 26.6 Gyr -1 (equivalent to a period of ~ 0.24 Gyr). Thus using Eq. ( fil| ) 



F (t) ~0.9JV, 



fa(x)\ 


2 




I kpc ; 


km s- 1 


V 10 Gyr) 



sat 

If now we assume that the progenitor systems are similar to present-day dwarf ellipticals, then using Eq. ([46j) we find for the 

-"3 

t \ A \ 5.1 x 10 2 , 100 x 10 r L sat, 



F&{t) ~ {TOGfr ) X \3.0x 10< 10 x 10 s L sat! (48) 
For t ~ 10 Gyr, the number of streams expected in the solar neighbourhood is therefore in the range 

F e ~ 300 - 500. (49) 



Fuchs & Jahreifi (1998) have obtained a lower limit for the local mass density of spheroid dwarfs of 1 x 10 s M kpc -3 . 
We may use this estimate to derive the mass content in subdwarfs of an individual stream in a volume of 1 kpc 3 centered on 
the Sun: 

F M (t)~ M ^^°J™ lkpcS l (50) 
f (t) 

Thus with our previous estimate for the filling factor 

/10Gyr\ 3 / 1.9xlO 2 M , for TO 7 L sat, 
Hm ~ \~^) X \ 3.3 x 10 2 M , for 10 8 L sat. (51) 

Therefore, after 10 Gyr, each stream contains Fm ~ (200 — 350) M in subdwarf stars, depending on the orbital parameters 
of the progenitors and their initial masses. 

Since the halo stars in the solar neighbourhood have one-dimensional dispersions <r b s (^) ~ 100 — 150 kms , in order to 
distinguish kinematically whether their distribution is really the superposition of ~ 300 — 500 individual streams of velocity 
dispersion <7 st (i>) we might require that 

(T st {v) < — — , (52) 
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where the factor 1/27 would ensure a ~ 3<r distinction between streams. Using our previous estimate of Fq this condition 
becomes <r B t(v) < cr obs (u)/(20 — 24), and thus a st (v) < 5 kms . Currently the observational errors in the measured velocities 
of halo stars are of order 20 kms -1 , and thus there is little hope to distinguish at the present day all the individual streams 
which may make up the stellar halo of our Galaxy. Since intrinsic velocity dispersions for streams originating from 10 7 — 1O 8 L0 
objects are of the order of 3 — 5 kms -1 after 10 Gyr, it should be possible to distinguish such streams with the astrometric 
missions SIM and GAIA, if they reach their planned accuracy of a few kms -1 . Even with an accuracy of 15 kms -1 per velocity 
component, streams are predicted to be marginally separated. The clumpy nature of the distribution should thus be easily 
distinguishable in samples of a few thousand stars. One way of identifying streams which are debris from the same original 
object, is through clustering in action or integrals of motion space (Helmi, Zhao & de Zeeuw 1998). 



5 AN OBSERVATIONAL APPLICATION 

Majewski et al. (1994) discovered a clump of nine halo stars in a proper motion survey near the NGP (Majewski 1992), which 
appeared separated from the main distribution of stars in the field. They measured proper motions, photometric parallaxes, 
F magnitudes and (J — F) colours for all nine stars and radial velocities for six of them. For these six stars we find for the 
mean velocity v v — —152 ± 23 kms -1 , vr = —260 ± 18kms -1 and v z — — 76 ± 18kms -1 , and for the velocity dispersions 
(j{v v ) = 99 ± 33 kms -1 , o{vr) — 100 ± 24kms -1 and o(v z ) — 35 ± 24kms -1 . If the dispersions are computed along the 
principal axes, we find <r(vi) = 29 ± 20 kms -1 , a(v2) — 68 ± 94 kms -1 , a(vg) — 125 ± 5 kms -1 . 

Since the mean velocities are significantly different from zero, the group of stars can not be close to any turning point of 
their orbit. The only way to understand the large observed dispersions, in particular of a{v3), if the stars come from a single 
disrupted satellite, is for the group to consist of more than one stream of stars. We believe that this may actually be the case. 
By computing the angular momenta of the stars we find they cluster into two clearly distinguishable subgroups: l}^ — —784 
and a^\L z ) — 299, and L z = —2180 and a^ 2 \L z ) = 313 in kpc kms -1 . If we accept the existence of two streams as a 
premise, we may compute the velocity dispersions in each of them. We find for the stream with 4 stars 

cr (1) («i) = 25 ± 25, a w {v 2 ) = 43 ± 62, a (1) (v 3 ) = 100 ± 45, 

while for the stream with 2 stars 

cr (2) («i) = 3 ±4, g (2) {v 2 ) = 25 ± 21, a {2) (v 3 ) = 89 ± 64, 

all in kms -1 . These results are consistent at a 2a level with very small 3-D velocity dispersions, as expected, if indeed these 
are streams from a disrupted satellite. 

With this interpretation of the kinematics of this group, we can estimate the mass of the progenitor and its initial size 
and velocity dispersion. Galaxies today obey scaling laws of the Faber-Jackson or Tully-Fisher type. If we assume that the 
original satellite was similar to present-day dwarf ellipticals, then we may use Eq. ( [IB] ) to derive a relation between the initial 
dispersion in the z-component of the angular momentum and initial velocity dispersion of the progenitor 



af(L z ) = af(v)Rt po + 0.0375 2 #£-<r? '», (53) 

where i? apo is the apocentric distance of its orbit. Under the assumption that L z is conserved, we can derive (Ji(v) by replacing 
in the previous equation the observed values of L z , a(L z ) and an estimate of -R apo . We obtain the latter by orbit integration in 
a Galaxy model, which includes a disk, bulge and halo and find i? apo ~ 12 kpc. Our estimate for the initial velocity dispersion 
of the progenitor is then 

cn{v) ~ 48 kms -1 , (54) 
which in Eqs. ( fl"!] ) and (JIB]) yields for its initial luminosity and size 

L~2xlO 8 L , i?~lkpc. (55) 

We estimate that the relative error-bars in these quantities are of order 50%, if measurement errors and a 50% uncertainty in 
the apocentric distance are included. 

In summary, if indeed these stars come from a single disrupted object, we must accept that the first six stars that were 
detected (Majewski et al. 1994) are part of at least two independent streams. This seems reasonable, since two streams can 
be indeed be distinguished, and the velocity dispersions, in each stream are very small. Moreover, a disrupted object with 
the properties just derived (luminosity, initial size and velocity dispersion), would fill its available volume rapidly, producing 
a large number of streams. In view of our explanation, a number of stars from the same disrupted object but with positive 
z-velocities should also be present in the same region, since phase-mixing allows streams to be observed with opposite motion 
in the R and/or z directions. Candidates for such additional debris should have similar v v , since L z is conserved during 
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phase- mixing. By simple inspection of Figure 1(a) in Majewski et al. (1994), other stars can be indeed found, with similar v v 
but opposite vr and v z . 



6 DISCUSSION AND CONCLUSIONS 

We have studied the disruption of satellite galaxies in a disk + halo potential and characterised the signatures left by such 
events in a galaxy like our own. We developed an analytic description based on Liouville's theorem and on the very simple 
evolution of the system in action-angle variables. This is applicable to any accretion event if self-gravity is not very important 
and as long as the overall potential is static or only adiabatically changing. Satellites with masses up to several times 10 9 Mq 
are likely to satisfy this adiabatic condition if the mass of the Galaxy is larger than several times 10 10 Mq at the time of infall 
and if there are no other strong perturbations. Even though have not studied how the system gets to its starting point, it seems 
quite plausible that in this regime dynamical friction will bring the satellites to the inner regions of the Galaxy in a few Gyr, 
where they will be disrupted very rapidly. Their orbital properties may be similar to those found in CDM simulations of the 
infall structure onto clusters, where objects are mostly on fairly radial orbits (Tormen, Diaferio & Syer 1998); this is consistent 
with the dynamics of solar neighbourhood halo stars. Their masses range from the low values estimated observationally for 
dwarf spheroidals to the much larger values expected for the building blocks in hierarchical theories of galaxy formation. 

We summarize our conclusions as follows. After 10 Gyr we find no strong correlations in the spatial distribution of a 
satellite's stars, since for orbits relevant to the bulk of the stellar halo this is sufficient time for the stars to fill most of their 
available configuration volume. This is consistent with the fact that no stream-like density structures have so far been observed 
in the solar neighbourhood. On the contrary, strong correlations are present in velocity space. The conservation of phase-space 
density results in velocity dispersions at each point along a stream that decrease as 1/t. On top of the secular behaviour, 
periodic oscillations are also expected: at the turning points of the orbit the velocity dispersions, and thus the mean density of 
the stream, can be considerably enhanced. Some applications of this density enhancement deserve further study. For example, 
the present properties of the Sagittarius dwarf galaxy seem difficult to explain, since numerical simulations show that it could 
have been disrupted very rapidly given its current orbit (Johnston, Spergel & Hernquist 1995; Velazquez & White 1995). This 
puzzle has led to some unconventional suggestions to explain its survival, like a massive and dense dark matter halo (Ibata & 
Lewis 1998) or a recent collision with the Magellanic Clouds (Zhao 1998). However, since the densest part of Sagittarius seems 
to be near its pericentre, it could be located sufficiently close to a 'caustic' to be interpreted as a transient enhancement. 
Sagittarius could simply be a galaxy disrupted several Gyr ago (c.f. Kroupa 1997). 

If the whole stellar halo of our Galaxy was built by merging of N sat similar smaller systems of characteristic size <j(x) 
and velocity dispersion cr(v), then after 10 billion years we expect the stellar distribution in the solar neighbourhood to be 
made up of Fq streams, where 

F ~O.9AU f^y^L.. 

V kpc j km s 1 

For satellites which obey the same scaling relations as the dwarf elliptical galaxies, this means 300 to 500 streams. Individually, 
these streams should have extremely small velocity dispersions, and inside a 1 kpc 3 volume centered on the Sun each should 
contain a few hundred stars. Since the local halo velocity ellipsoid has dispersions of the order of 100 kms -1 , 3-D velocities 
with errors smaller than 5 kms -1 are needed to separate unambiguously the individual streams. This is better by a factor 
of four than most current measurements, which would, however, be good enough to give a clear detection of the expected 
dumpiness in samples of a few thousand stars. The combination of a strongly mixed population with relatively large velocity 
errors yields an apparently smooth and Gaussian distribution in smaller samples. Since the intrinsic dispersion for a stream 
from an LMC-type progenitor is of the order of 3 — 5 kms -1 after a Hubble time, one should aim for velocity uncertainties 
below 3 kms -1 . With the next generation of astrometric satellites, (in particular GAIA, e.g. Gilmore et al. 1998) we should 
be able to distinguish almost all streams in the solar neighbourhood originating from disrupted satellites. 

Our analytic approach is based on Liouville's Theorem and the very simple evolution of the system in action-angle 
variables. Although the latter is likely to fail in the full merging regime, the conservation of local phase-space density will 
still hold. It will be interesting to see how this conservation law influences the final phase-space distribution in the merger of 
more massive disk- like systems. These are plausible progenitors for the bulge of our Galaxy in hierarchical models. 
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APPENDIX A: SPHERICAL POTENTIAL 
Al 2-D case 

For a spherical potential "l>(r), the Hamiltonian is separable in spherical coordinates and depends on the actions and Je 
only through the combination J^ + Jg = L. We therefore may choose a system of coordinates which coincides with the plane of 
motion of the system, reducing the problem to 2-D. The position of a particle is given by its angular ip and radial r coordinates 
on that plane. In that case, we have 

m 

(Al) 



L = J i ,=p i ,, J r = - [ dr-J2(E-${r))r 2 -L 2 
7r / r 



where L is the total angular momentum of the particle, E its energy and r\ and r-2 are the turning points in the radial direction 
of motion. The action J r cannot be computed analytically in general for an arbitrary potential. However, Eq. (Al) defines 
an implicit function g — g(E,L, J r ) = 0, which we can differentiate to find the frequencies of motion and their derivatives. 
These are needed to compute the elements of the matrix ®(t) (Eq. (^)) and to obtain the time evolution of the distribution 
function. 

To simplify the computations, we assume that the variance matrix |^] in action-angle variables is diagonal at t = 0: 
G wij = (JuSij. The evolution of the system in phase space is obtained through the product 0(t) V°„0(t), which yields the 
following variance matrix a w (t)ij = {i, j} at time t 



>{t) = 



(Til —(711^x1^ — 0"11^12^ 

{1,2} (722 — C"22^'l2i — 022^22^ 

{1,3} {2,3} ffntfuV+osafiia * + 



{1,4} {2,4} 



{3,4} 



^22^22 ^ ~t~ ^"44 



/ 2 ,2 



in action-angle variables, with fiy = dQi/dJj. Subindices {1} and {3} refer to directions associated to ip, such as for example 
(j>^ and J^, whereas {2} and {4} are related to r. 

We find the properties of the debris in configuration and momenta space by transforming the action-angle coordinates 
locally around x with the matrix T _1 . Its elements are the second derivatives of the characteristic function W(q, J): 



WjjJ p 

J r, 



(A2) 



with J a — — W„ } W 00 and J„ = W ? j , and has the following form for a spherical potential in 2-D 



qj 

tl2 
t>-2 



ti-2 



tl3 

1 

t-4?, 



tl4 

t-24 


t-44 



(A3) 



with 

tl2 = 
^22 = 

t42 - 

and 

h(r) = 



+ 



2 



h(r) d 2 W 



fir dLdJ, 
h(r) 2 W + fi_ 

^iif dJ r pT* 

h(r) 



-$'(r) + — , p r = 





d 2 W d 2 W 




tl3 = 


dL 2 1 dLdJr ' 


ti4 — 




d 2 w d 2 w 




t23 = 


OLdJr ' dJ? t43 ' 


t24 = 








t43 = 


-a: 


t44 — 



d 2 W 



OLdJr fir 
2 W p r 

Pr 



2(E-*(r))-^, 



t As we mentioned in Section 4.2, a is in fact the inverse of the covariance matrix. However we refer to a as the variance matrix. 
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where all functions are evaluated at x. Therefore the variance matrix in (x, p) is 
a u = (0(f)T- 1 ) t a° m (0(f)T- 1 ), 
so that, by substituting 



(A4) 





(Til 


(TllA 




(TiiB anC 




{1,2} onA 2 - 


r a 22 D 2 


+ (744 £42 


AauB + D022E + f420"44i43 AanC + D022F + t^iitii 




{1,3} 


{2,3} 




OllB 1 + (J22E 2 + CT33 + CT44^ 3 BtJuC + E022F + 1 43 O 44t 44 




{1,4} 


{2,4} 




{3, 4} (JnC 2 + G22F 2 + (744*44 


and where 








A = 


= tl2 — fl'i2t42t, 


B = 


tl3 — fl'nt 


— r2' 12 *43i, C = tl4 — Q'i2t44t, 


D = 


- t22 — 5^22^42^, 


E = 


t23 — fi'l2^ 


— ^22^43^, F = t24 — fl^tni. 



In general, one is more interested in the characteristics of the debris in velocity space, rather than in momenta space. Thus 
we transform the variance matrix according to <7 W = Tl^, v cr w T-p^ v , with 



10 

10 

ity r 

1 



(A5) 



The diagonalization of the variance matrix yields the values of the dispersions along the principal axes and their 
orientation: two of its eigenvalues increase with time, whereas the other two decrease with time. To understand the directly 
observable properties of the debris we examine what happens around a particular point x(t) in configuration space located 
on the mean orbit of the system. This is equivalent to studying the velocity submatrix 



E 2 + (733 + (744^3) r(BanC + E022F + tuaatii) 



{1,2} 



CllC +(7 2 2-F +(744^44 



For example, by diagonalising the matrix cr ro (v) we obtain the directions of the principal axes of the velocity ellipsoid at 
the point x(t), and their dispersions. Its eigenvalues are the roots of det[<7 ro (u) — AX] = 0. For t 3> t OI h 



A * = 2 



a n fi'u — — — ( — - 



L 



(fir) 



(7llfi'l2 +(722^22 



fi'l2 ^ ( ^1/' o 



+ 



± v^R 



(A6) 



for i — 1,2, and where 



/ / 22/ J " 
+ (722 | f2 12 ( 7 



+ 4r 

Therefore 
A1A2 = 



o' 2 1 n' 2 
C11SJ12 +(722"22 



(7iif2'i 2 ( fi'u ^y 2 - ( 0^ — — 



+ (722^22 f ^12 ~~ ^2 



(A7) 



<7ll<72 2 r , 



Al + A 2 = t 2 r 2 \(7u(Cl'n + Q.^Azf + (722(^12 + ^M^) 2 ] + t 2 



(711^1 



- (72 2 fi 2 ; 



Since cr(t)j) = y^l/Xi both velocity dispersions decrease on the average as 1/t. The principal axes of the ellipsoid rotate as 
time passes by, not being coincident with any particular direction. 



A2 3-D treatment 

As we discussed in Section (4.2), the problem of the disruption of the system and its evolution in phase-space is really a 3-D 
problem, since our initial satellite had a finite width in all directions. Since we just discussed in great detail what happens in 
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the 2-D case and the way of proceeding once more dimensions are added is the same, we will simply outline our main results, 
focusing on what happens to the velocity submatrix. 

If we assume that the system had initially a diagonal variance matrix in action-angle variables, the velocity submatrix at 
time t is 



a m (v) = T^„S(t) T p _„, 
with 

s(t) = (Wx/W-, 1 -m'w-j)^°jwjjwzj - tn'wrj) + w-jMwr, 1 , 



(A8) 



(A9) 



where Wjj is the matrix whose elements are the second derivatives of the characteristic function with respect to the 
actions, W q j the matrix that contains the second derivatives of W with respect to the coordinates q and the actions J, 
and fiy = dQi/dJj. Note that, since the potential is spherical fi' and Wjj have two equal rows. The initial variance matrix 
in action-angle space 

4 



We can compute the density at a later time at the point x(t) located on the mean orbit of the system by integrating 
/(x,v,t) = /oexp r ' 



-At(t)a ro (t)A ro (t) 



with respect to the velocities using the submatrix cr OT (u) 



p(x,t) = 



dv v dve dv r /(x, v, t). 



In the principal axes frame 
p(x,t) = f Q (2ir) 3/2 a vl (t)a V2 (t)a V3 (t)Eri 



ai 



V2a vi (t) 



Erf 



"2 



Erf 



(>■>, 



V2<r V3 (t) 



(A10) 



V2(T V2 (t)_ 

with ai,a,2,a3 the boundaries of the integration volume. For t 2> £ rb the error function tends to 1, and therefore 

p(x,t) = (2n f /2 f a vl (t)a V2 (t)a V3 (t), (All) 

that is equivalent to 

(A12) 



p(x,t) = (2TT) 3/2 f /y/\i\ 2 \3, 
where the A's are the eigenvalues of a^(v). With simple algebra it can be shown that 
A1A2A3 = det avj(v), 



which is readily computable from Eqs. (A8) and ( A9) 

det a^(v) = (det T p ^„) 2 (det W^ 1 ) 2 det[(Wjj - t n')V$(Wjj - t n') + <t% 
where 

det Tp^ v = r 2 sin 8 
and 



detW-, 1 



Pr Ve 

L 



(A13) 



(A14) 



(A15) 



(A16) 



The remaining determinant in Eq. (A14) for t 2> iorb is 

C"33(o"ll + 022X0-44 + Cr55)(Sl' li r2 33 — Q[ 3 2 ) 2 t 4 , 

so that finally 

(27t) 3/2 / Q r L 



p(x,t) = 



\Q' U Q,' 33 - 0' 13 2 I % /<T33(<Tll+a 2 2)(<T44+a 55 ) r 2 sin 61 |W« I t2 ' 

Let us recall that an = for i = 1..3 and an — l/a 2 . for i = j + 3 = 4. .6. 
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APPENDIX B: AXISYMMETRIC EDDINGTON POTENTIAL 

To exemplify and understand how the rupture of the spherical symmetry affects the characteristic scales of the system, we take 
a very simple Eddington potential $(r, 9) = <J>i(r) + 77 (/3 cos 8)/r 2 (Lynden-Bell 1962; 1994) which is separable in spherical 
coordinates. The third integral for this class of potentials is I3 = 5L 2 + ri((3cos8). The actions are computed from: 



Je = 



2tt 



de X 2(h- v (e)) 



J 2 . 



sin 2 6 7 



*W2(£-$ lW )_5. 



(Bl) 
(B2) 

(B3) 



The procedure outlined in Section 4.1 and Appendix A can also be applied to a system moving in this type of potentials. 
In particular we are interested in the behaviour of the density. By virtue of the previous discussion we only need to find the 



determinant of the variance matrix as in Eq. (A14), for this potential. Since Eqs. (A15) and (A.16) remain unchanged, we 
only focus on det[(WjjW~j - tSl'Wlj^ cr^ (WjjW"j - t(l"W~j) + a% For t > t orh the term with fl' will dominate with 

respect to Wjj, and the product t 2 (ft''W~j)^a^ > fl'W~j will dominate over cr°jQ Therefore 
det cr w (u) = (detTp_»„) 2 (detW^ J 1 ) 2 (detn't) 2 deta^, 
and so the density at the point x at time t is 

(27T) 3 / 2 /o 1 8h Q r 1 



(B4) 



p(x,t) 



(B5) 



ydet 7 ^ |detO'| dJe r 2 sin 9\p rP g\ t 3 
This expression is valid for a satellite described initially by a Gaussian distribution. The variance matrix at t = t° may 



be 



(i) diagonal in action-angle variables: 
deto-0 = 1/((70 1 <T0 2 (T^ 3 ) 2 , 

(ii) diagonal in configuration- velocity space: 



det a. 



PePr 1 I I 



{l 2 (dh/dJe) 2 a 2 a 2 



W 2 Wk.J / COS 2 ' 



sin 2 6 



W, 



;. + ( 



I'- 
ve cost 



r sin ( 



Wee \ 



where all functions are evaluated at (x , v ), and 



Wee 



and 



he 
pe ' 



he = -v'(0) + J, 



2 cos ( 



sin d 6 



Wrr = 



Pr 



2/a 



The expression for the determinant of the angle submatrix at t = may be simplified if the satellite is initially close to a 
turning point of the orbit. In this case the term WegW 2 r will be dominant and 

2 



deter rf 



he(e°)h r (r°) (dhy 1 



1 



2 

o v o v 



Note that the main differences with the spherical case are 

• the time dependence: t 3 instead of t 2 because of the increase in the dimensionality of the problem; 

• the dependence on the derivatives of the basic frequencies of motion: the same functional dependence detfi', but now 
with three independent frequencies and derivatives; 

• the inclusion of the term dls/dJe, which for the spherical case is simply L; 

• the form of pe = \l 1(1$ — 77(6*)) — J 2 / sin 2 6, which also includes the angular dependence of the potential. 



S This does not hold for the spherical case because det[(f2'W" J l )1' f TOO'W~ 7 1 ] oc det ST = 
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APPENDIX C: AXISYMMETRIC STACKEL POTENTIAL 



In this section we collect some basic properties of Stackel potentials and derive the density behaviour as a function of time, 
as in previous sections, from Liouville's Theorem and the evolution of the system in action-angle variables. Further details on 
Stackel potentials can be found in de Zeeuw (1985). 

Let us first introduce spheroidal coordinates (A, u, <p), where p is the azimuthal angle in the usual cylindrical coordinates 
(R, z, ip), and A and v are the two roots for r of 



R 



+ 



= 1, 



t — or r — 

where c 2 < v < a 2 < A. A potential is of Stackel form if it can be expressed as 
(\-c 2 )G(\)-(u~c 2 )G(v) 



V = -- 



where G(r) is an arbitrary function (r = A, v). In this case, the Hamiltonian becomes 



where the functions P and Q are 

„a A — v „ 5 



v — \ 



4(A-a 2 )(A-c 2 )' 



i(v — a?)(v — c 2 ) ' 



(CI) 



(02) 



(03) 



(04) 



Three isolating integrals of motion can be found (E, I2, I3), and the system is separable since the equations of motion can be 
written as 



1 



G(r) 



(05) 



1T 2(r-a 2 ) 
and 

p v = L z = V2h- (06) 

To represent the Galaxy we may choose a superposition of two Stackel potentials: a disk plus a halo component 

V = kV disk + (l-k)V ha i , (07) 

where k represents the mass fraction of the disk with respect to the total mass of the Galaxy. Since the coordinates used 
for the halo and the disk have to be the same, this introduces a relation between the characteristic parameters (eid,Cd) and 
(ah, Ch) of the Stackel potentials. It can be shown that the potential 



V(X,u,q) 



-GM 



+ 



VX + ^fv ^/\ + q + yjv + q 



(08) 



where g is related to the flattening of the halo component, provides a good description yielding a flat rotation curve with 
similar properties to that of our Galaxy (Batsleer & Dejonghe 1994). The function G(t) in Eq. (C2) is 



G(r) = GM 



k 



+ 



^/t + c y/r + q + c 



(09) 



For the characteristic parameters we choose ad = 2, Cd = 1, ah/ch = 1.01 (giving a rather spherical halo), k = 0.12 and 
M = 5 x 10 11 M Q . 

In order to obtain the evolution of the mean density of debris as a function of time in a Stackel potential we use the 
results of Section 4 and of Appendix A and B. From Eqs. (A12) and (A13) the density is proportional to the determinant 
of the velocity submatrix. Since the Hamiltonian is separable in spheroidal coordinates, to obtain the density in cylindrical 
(or spherical) coordinates we need to multiply Eq. (B4) by the determinant of the matrix that performs the transformation 
between the two sets of coordinates. Thus 



det (v) — [det T p ^„ det T PT ^ Pcyl det W ? j det tt'tj 2 det cr^, 
where 

det T p _„ = R, det T Pr ^ Pcyl det W^ 1 = ^ ~ ^ VxV " 



o 9/3 O dh 

OJx OJ v 



(CIO) 

(on) 



The mean density at time t at the point x on the mean orbit of the system becomes 
© RAS, MNRAS 000, 



28 A. Helmi and S.D.M. White 
{^f ,2 h i 



p(x,t) 



OJx oJ v 



(C12) 



^deta° \ Vx v u \\\-v\R |detfi'| t 3 
This expression is valid for a satellite described initially by a Gaussian distribution. The variance matrix at t — t° may be 

(i) diagonal in action-angle variables: 
det<T0 = l/(a r j >1 a l p 2 a l f, 3 ) 2 , 

(ii) diagonal in configuration-velocity space. If the satellite is initially close to a turning point of the orbit then 

1 2 

h\h v A — v 1 



det ctj 



oJx oJ v 



P 3 Q 3 o^gI 



(C13) 



where all functions are evaluated at (x , v ), and 

u r, dPr , 
h T — 2p T ——, t = X,v. 
or 
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